基于MODIS数据的青藏高原1-km分辨率
逐日地表温度数据集(2000–2020)
徐迅澎1, 2, 3，张 煜1, 2, 3，张玉常4，计璐艳1, 2，唐海蓉1, 2, 3*
1. 中国科学院空天信息创新研究院，北京
100094；
2. 中国科学院空间信息处理与应用系统技术重点实验室，北京
100190；
3. 中国科学院大学电子电气与通信工程学院，北京
101408；
4.
火箭军某军代室，北京 100000
摘 要：遥感数据在空间、时间上具有较强相关性和连续性，因此遥感图像时间序列具有低秩属性。在本文中，我们考虑利用低秩张量补全的方法修复图像。首先，我们对MODIS的地表温度数据进行预处理，通过时空插值的方式将被云遮挡导致的缺失值进行初步补全。接着，我们将地表温度时间序列数据视作一个三阶时空张量，在时空张量的时间维度上引入傅里叶变换，将其转换为空间-频率张量。对此张量进行奇异值分解、高斯低通滤波后再通过逆傅里叶变换为空间-时间张量。最后使用交替方向乘子法迭代运算，对缺失张量进一步优化。对于该方法的精度，我们采用模拟实验的方法进行验证，人工添加掩膜并进行恢复，恢复结果的平均绝对误差（Mean Average Error）在2.1–4.9 K内。本数据集包括青藏高原2000–2020年每日的以下数据：（1）对MOD11A1、MYD11A1产品中被云遮挡区域进行优化的地表温度数据（MOD11A1_QTP_PART，MYD11A1_QTP_PART）；（2）被云遮挡区域优化后的MOD11A1、MYD11A1产品，即优化后的地表温度数据（MOD11A1_QTP_Temp，MYD11A1_QTP_Temp）；（3）原始的MOD11A1、MYD11A1产品（MOD11A1_QTP_ORIGIN，MOD11A1_QTP_ORIGIN）。所有数据空间分辨率为1 km，以整数型数据格式进行存储，像元值表示地表的热力学温度，比例系数为0.02，单位为K。数据集存储为.tif格式，可以通过ENVI、ArcGIS等遥感软件直接打开和处理。
关键词：青藏高原；逐日地表温度；1 km；2000–2020；MODIS
DOI: https://doi.org/10.3974/geodp.2023.03.03
CSTR: https://cstr.escience.org.cn/CSTR:20146.14.2023.03.03
数据可用性声明：
本文关联实体数据集已在《全球变化数据仓储电子杂志（中英文）》出版，可获取：
https://doi.org/10.3974/geodb.2023.10.02.V1或https://cstr.escience.org.cn/CSTR:20146.11.2023.10.02.V1.
1 前言
地表温度（Land Surface
Temperature，LST）是指地表在近红外和热红外波段反射和辐射的能量，是用于描述地表热态的重要指标之一，该指标在气候变化、生态环境和农业生产等领域具有重要的研究价值[1–3]。使用遥感数据来获取地表温度能够快速、全面、准确地获取大范围内的地表温度分布情况，从而能够有效地指导相关领域的决策和规划。并且随着遥感技术的不断发展和卫星数据的不断更新，地表温度数据的分辨率和精度得到了大幅提升，有越来越多的地表温度产品被生产出来，供研究人员们使用。
青藏高原地区是我国重要的区域之一，其广阔的面积、丰富的资源、独特的地理环境以及特殊的气候条件，使其一直是地理、气候和环境科学领域中的研究热点之一。近些年来，随着全球气候变化加剧，青藏高原地表温度异常变化问题引起了广泛关注和研究。因此，对青藏高原地表温度进行研究，关系到气候变化、水资源管理、生态环境保护等多方面的问题，具有重要的理论与实践意义。
然而，由于该地区的地形特殊、气候环境复杂，给该地区的遥感数据获取和处理带来了较大的挑战。其中一个挑战就是云遮挡的问题，当被云遮挡的区域过多时，可能会对地表温度数据的精度和可用性造成较大的影响，影响其应用价值。因此，如何提高地表温度数据的质量和可用性是当前亟待解决的核心问题之一。
低秩张量补全方法是一种有效的技术，可以用于遥感图像中的云修复[4–7]。该方法能够充分利用数据的低秩属性，从不完整的观测中进行张量补全，从而恢复遥感数据中被云层遮挡的信息。该方法能够有效地保留细节信息，提高数据的质量和应用性。
本文采用MOD11A1 V6和MYD11A1 V6产品的逐日地表温度数据和时空联合低秩张量补全的方法[8]，对青藏高原地区逐日地表温度数据进行缺失的数据补全并裁剪[9]。最终生产出2000–2022年间青藏高原地区的逐日无云地表温度数据集。该数据集具有广泛的研究和应用价值，能够促进青藏高原地区的气候研究、生态环境评估等相关领域的发展。
2 数据集元数据简介
《基于MODIS数据的青藏高原1-km分辨率逐日地表温度优化数据集（2000–2020）》[10]的名称、作者、地理区域、数据年代、时间分辨率、空间分辨率、数据集组成、数据出版与共享服务平台、数据共享政策等信息见表1。
3 数据研发方法
3.1 数据源
本数据集采用的数据源来自MOD11A1 V6[1]和MYD11A1 V6[2]，它们是基于MODIS的逐日地表温度产品，能够提供陆地表面的日间和夜间地表温度数据，该数据产品覆盖了全球，包括青藏高原地区。其中，MOD11A1来自于Terra卫星，而MYD11A1来自于Aqua卫星，两颗卫星过境的时间不同，可以分别获取当地上午和下午的遥感数据。
表1 《基于MODIS数据的青藏高原1-km分辨率逐日地表温度优化数据集（2000–2020）》元数据简表
|
条 目 |
描 述 |
|
|
数据集名称 |
基于MODIS数据的青藏高原1-km分辨率逐日地表温度优化数据集（2000–2020） |
|
|
数据集短名 |
MODIS_QTP_Temp |
|
|
作者信息 |
徐迅澎，中国科学院空天信息创新研究院，中国科学院空间信息处理与应用系统技术重点实验室，中国科学院大学电子电气与工程学院，xuxunpeng21@mails.ucas.ac.cn 张煜，中国科学院空天信息创新研究院，中国科学院空间信息处理与应用系统技术重点实验室，中国科学院大学电子电气与工程学院，zhangyu217@mails.ucas.ac.cn
计璐艳，中国科学院空天信息创新研究院，中国科学院空间信息处理与应用系统技术重点实验室，jily@mail.ustc.edu.cn 唐海蓉，中国科学院空天信息创新研究院，中国科学院空间信息处理与应用系统技术重点实验室，中国科学院大学电子电气与工程学院，tanghr@aircas.ac.cn |
|
|
地理区域 |
青藏高原 |
|
|
数据年代 |
2000–2020年 |
|
|
时间分辨率 |
1天 |
|
|
空间分辨率 |
1 km |
|
|
数据格式 |
.tif |
|
|
数据量 |
138 GB（压缩后） |
|
|
数据集组成 |
包括青藏高原2000–2020年逐日的以下数据：（1）对MOD11A1、MYD11A1产品中被云遮挡区域进行优化的地表温度数据（MOD11A1_QTP_PART，MYD11A1_QTP_PART）；（2）被云遮挡区域优化后的MOD11A1、MYD11A1产品，即优化后的地表温度数据（MOD11A1_QTP_Temp，MYD11A1_QTP_Temp）；（3）原始的MOD11A1、MYD11A1产品（MOD11A1_QTP_ORIGIN，MOD11A1_QTP_ORIGIN），每个目录下数据的命名规则为YYYYDDD.tif，其中YYYY代表年份，DDD代表某年的第多少天，例如2020001.tif |
|
|
基金项目 |
中华人民共和国科学技术部（2019QZKK0206，31400） |
|
|
数据计算环境 |
Python |
|
|
出版与共享服务平台 |
全球变化科学研究数据出版系统
http://www.geodoi.ac.cn |
|
|
地址 |
北京市朝阳区大屯路甲11号100101，中国科学院地理科学与资源研究所 |
|
|
数据共享政策 |
（1）“数据”以最便利的方式通过互联网系统免费向全社会开放，用户免费浏览、免费下载；（2）最终用户使用“数据”需要按照引用格式在参考文献或适当的位置标注数据来源；（3）增值服务用户或以任何形式散发和传播（包括通过计算机服务器）“数据”的用户需要与《全球变化数据学报（中英文）》编辑部签署书面协议，获得许可；（4）摘取“数据”中的部分记录创作新数据的作者需要遵循10%引用原则，即从本数据集中摘取的数据记录少于新数据集总记录量的10%，同时需要对摘取的数据记录标注数据来源[11] |
|
数据和论文检索系统 |
DOI，CSTR，Crossref，DCI，CSCD，CNKI，SciEngine，WDS/ISC，GEOSS |
在本文中，我们选取其中的LST_Day_1km波段进行恢复，该波段是MOD11A1 V6和MYD11A1 V6产品的一个重要波段，用于获取白天地表温度信息，空间分辨率为1 km，空间分辨率为1天，一景数据的尺寸为1200×1200。LST_Day_1km波段的数据有效范围为7,500–65,535，比例系数为0.02，无效值则被设为0。
3.2 数据预处理
青藏高原多数地区常年处在多云的状态，这会导致地表温度数据中存在许多常年被云遮挡的噪音。如果直接使用联合时间维FFT的低秩张量补全算法，会将常年出现的云视作低频成分，从而代替本应复原的地物位置，最终导致修复结果中出现许多黑色条带，影响数据的可用性。
因此，我们首先需要对地表温度数据进行预处理，通过时空插值的方式将大部分缺失的有云位置补全，示意图如图1所示。首先，先将多时相的地表温度数据按时间维进行排列，记作
，其中
和
是数据的空间维大小，
是时间序列长度。其次，将地表温度时间序列
划分为多个大小为
的小窗口，对各个窗口内的有效值进行平均处理，将其作为窗口下采样后的值，得到下采样后的地表温度时间序列
。通过最小二乘法，逐像素求解出原始地表温度时间序列与降采样后的地表温度时间序列有效值之间的线性关系。最后，利用求解得到的系数与下采样后的地表温度时间序列，对有云遮挡的地方的数据恢复，最终得到预处理后的地表温度时间序列
。经过数据预处理后，我们能够把数据中的大部分云给去除。此时，云出现的频率大大减少，便于我们进行后续的张量补全操作。
图1 预处理操作示意图
由于研究所用的数据时间序列较长，利用无云期间的地表温度值进行时空插值的过程可以提高图像获得的信息的稳定性和精度。此外，以100×100的窗口为单位进行降采样，也能够保证数据的局部空间一致性，避免出现突变的情况。同时，使用最小二乘法来建立原始地表温度时间序列与降采样后的地表温度时间序列有效值之间的线性关系能够更加准确地预测缺失数据点的值，避免了缺失点过多导致的误差累计的问题。
3.3 算法原理
由于遥感数据在空间、时间上具有较强的相关性和连续性，因此遥感图像时间序列，具有低秩属性，可以利用低秩张量补全的思想，达到修复图像的目的。在本文中，我们采用不同的分解方法来处理空间和时间维度，引入傅里叶变换来滤波时间维度，并根据时间维度的频谱属性自适应选择权重，应用到空间维度的低秩矩阵补全中。除此之外，我们还利用频域中的共轭对称性，加快了算法的运算速度。本文方法突出了时间维度上地物变化引起的低频成分，抑制由云引起的高频噪声部分，实现了时间空间联合低秩补全。
3.4 技术路线
《基于MODIS数据的青藏高原1-km分辨率逐日地表温度优化数据集（2000–2020）》研发的技术路线如图2所示。
3.4.1 时间维度FFT
傅里叶变换变换是将时域信号投影到一组正交三角函数基上的方法，适用于序列数据分解和处理。我们在张量的时间维度上引入傅里叶变换，将时域信息转到频域来处理，如下所示：
(1)
式中，
表示傅里叶变换后的张量，我们称为空间-频率张量。
3.4.2 时间维滤波
时间维度傅里叶变换后，时间序列地表温度数据的频谱可以分为低频成分和高频成分。低频成分对应于时间上缓慢变化或不变的部分，而高频成分对应于时间上变化明显的部分，比如云、噪声等。我们将在时间维方向上进行高斯滤波，保留其主要低频成分，削弱云、噪声等因素的影响。高斯滤波函数为
，其中
为像元个数，
为定义的常数。
3.4.3 空间维自适应低秩
为了达到时空和空间联合低秩，更好地选择和保存图像中的地物信息，对张量进行频域调制的自适应权重的低秩处理。对不同切片
根据矩阵信息的重要性，来决定每个切片低秩处理的阈值
：
(2)
(3)
3.4.4 时间维度iFFT
经过上面三步，我们再进行逆傅里叶变换，能够得到解的形式
。
将被云覆盖位置的值进行更新，直到更新值小于阈值后，我们就能够得到无云的LST时间序列。
图2 数据集制作技术路线图
4 数据结果与验证
4.1 数据集组成
本数据集按照卫星的不同分为6个目录：（1）对MOD11A1、MYD11A1产品中被云遮挡区域进行优化的地表温度数据（MOD11A1_QTP_PART，MYD11A1_QTP_PART）；（2）被云遮挡区域优化后的MOD11A1、MYD11A1产品，即优化后的地表温度数据（MOD11A1_QTP_Temp，MYD11A1_QTP_Temp）；（3）原始的MOD11A1、MYD11A1产品（MOD11A1_QTP_ORIGIN，MOD11A1_QTP_ORIGIN），每个目录下数据的命名规则为YYYYDDD.tif，其中YYYY代表年份，DDD代表某年的第多少天，例如2020001.tif.
4.2 数据结果
在本文中，我们选取MOD11A1_QTP_Temp产品中的2020001.tif的数据进行展示。图中黑色的部分为青藏高原之外的区域，在数据中用0表示。图片灰色部分为本文的研究区域，有效值范围为7,500–65,535。
我们选择在纳木错地区附近进行数据的定性分析。图4和图5展示了纳木错地区在2000–2009年和2010–2020年期间的地表温度数据的散点图。在图中，蓝色部分表示原始
图3 2020001.tif数据展示
图4 纳木错区域地表温度数据散点图（2000–2009）
图5 纳木错区域地表温度数据散点图（2010–2020）
产品中的有效值，我们将其保留。红色部分表示原始产品中的缺失数据，我们使用算法进行恢复。从下图可以看出，通过我们的算法恢复的数据基本上符合地表温度在一年内的变化趋势。
此外，我们还绘制了纳木错地区年平均地表温度的折线图，如图6所示。该折线图展示了纳木错地区在20年间的变化趋势。可以观察到，地表温度的变化相对平稳，相邻年份之间的变化约为1 K。
图6 纳木错区域2000–2020年年平均地表温度变化趋势图
4.3 数据结果验证
由于缺乏真实数据，我们使用模拟实验来验证数据集的恢复精度。我们以MOD11A1 V6产品2020年一整年的数据为例，随机挑选不同日期、不同位置的8块无云区域，为这些区域手动添加掩膜，使用0值代替原本产品中的地表温度信息，并用本文所提出的方法对确实的地表信息进行恢复。最终，我们使用评价指标对区域的恢复值进行评估，评估结果如表2所示。
表2 手工增加掩膜区域地表温度恢复效果表（无量纲）
|
参数 |
区域1 |
区域2 |
区域3 |
区域4 |
区域5 |
区域6 |
区域7 |
区域8 |
|
MAE |
3.013,4 |
3.812,5 |
4.912,9 |
4.333,4 |
4.716,9 |
2.806,4 |
2.112,0 |
3.624,1 |
|
RMSE |
3.992,6 |
4.553,3 |
6.164,6 |
5.402,0 |
5.659,0 |
3.582,7 |
2.689,3 |
4.731,6 |
|
R |
0.789,0 |
0.748,9 |
0.641,2 |
0.373,4 |
–0.810,3 |
0.554,4 |
0.694,2 |
0.408,1 |
除此之外，本文还将本产品与已有其他产品进行了比较，具体的比较结果如表3所示。其中，产品1是2020年青藏高原地区Landsat时序地表温度产品[12, 13]，产品2是中国区域1-km无缝地表温度数据集（2002–2020）[14–17]。
由于卫星过境的时间不同，使得各产品间无法直接进行比较。本文通过绘制相关地表温度的恢复情况变化趋势来比较本文方法的可靠性与准确性，变化趋势图如图7所示。
从上图可以看出，本文产品与其他产品在地表温度的恢复趋势上呈现出相同的趋势，在时间分辨率和空间精度上做到了较好的平衡。
表3 已有相关研究针对2020年纳木错区域地表温度恢复统计表（单位：K）
|
时间 |
01-01 |
01-17 |
02-02 |
01-31 |
11-16 |
12-02 |
12-18 |
|
产品1 |
240.00 |
279.30 |
239.90 |
308.70 |
292.80 |
283.50 |
278.50 |
|
产品2 |
276.80 |
280.32 |
278.92 |
304.64 |
285.36 |
280.08 |
279.66 |
|
本产品 |
268.96 |
275.84 |
278.24 |
296.16 |
284.32 |
276.64 |
274.56 |
图7 已有相关研究针对2020年纳木错区域地表温度恢复情况变化趋势
5 讨论和总结
在本文中，针对MOD11A1/MYD11A1
V6地表温度逐日产品的时间连续性，但存在大面积空间有效值缺失的特点，我们采用了时空插值以及时空联合低秩技术对缺失的有效值进行恢复。最终，成功地生产了青藏高原地区2000–2020年1-km逐日地表温度数据集。为了进一步验证产品的精度，我们与现有的其他两个产品进行了对比，结果表示我们的产品无论是时间分辨率还是空间精度上都有明显的优势，该数据集具有广泛的研究和应用价值，能够促进青藏高原地区的气候研究、生态环境评估等相关领域的发展。
作者分工：徐迅澎，张煜，计璐艳，唐海蓉对数据集的开发做了总体设计；张煜，计璐艳采集和处理了数据；徐迅澎，唐海蓉设计了模型和算法；张煜，计璐艳，张玉常做了数据验证；徐迅澎撰写了数据论文；张玉常提供了论文润色支持。
利益冲突声明：本研究不存在研究者以及与公开研究成果有关的利益冲突。
参考文献
[1] Wang, A. H., Zeng, X. B. Development of global hourly 0.5° land surface air temperature datasets [J]. Journal of Climate, 2013, 26(19): 7676–7691.
[2] Mostovoy, G. V., King, R. L., Reddy, K. R., et al. Statistical estimation of daily maximum and minimum air temperatures from MODIS LST data over the state of Mississippi [J]. GIScience & Remote Sensing, 2006, 43(1): 78–110.
[3] Xu, Y. M., Qin, Z. H., Shen, Y. Study on the estimation of near-surface air temperature from MODIS data by statistical methods [J]. International Journal of Remote Sensing, 2012, 33: 7629–7643.
[4] Ng, M. K-P., Yuan, Q. Q., Yan, L., et al. An adaptive weighted tensor completion method for the recovery of remote sensing images with missing data [J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(6): 3367–3381.
[5] Ji, T. Y., Yokoya, N., Zhu, X. X., et al. Nonlocal tensor completion for multitemporal remotely sensed images’ inpainting [J]. IEEE Transactions on Geoscience and Remote Sensing, 2018, 56(6): 3047–3061.
[6] Chen, Y., He, W., Yokoya, N., et al. Blind cloud and cloud shadow removal of multitemporal images based on total variation regularized low-rank sparsity decomposition [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2019, 157: 93–107.
[7] Lin, J., Huang, T. Z., Zhao, X. L., et al. Robust thick cloud removal for multitemporal remote sensing images using coupled tensor factorization [J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 1–16.
[8] Chen, Z. H., Zhang, P., Zhang, Y., et al. Thick cloud removal in multi-temporal remote sensing images via frequency spectrum-modulated tensor completion [J]. Remote Sensing, 2023, 15(5): 1230.
[9] 张镱锂, 李炳元, 刘林山等. 再论青藏高原范围[J]. 地理研究, 2021, 40(6): 1543–1553.
[10] 徐迅澎, 张煜, 计璐艳等. 基于MODIS数据的青藏高原1-km分辨率逐日地表温度优化数据集（2000–2020）[J/DB/OL]. 全球变化数据仓储电子杂志, 2023. https://doi.org/10.3974/geodb.2023.10.
02.V1. https://cstr.escience.org.cn/CSTR:20146.11.2023.10.02.V1.
[11] 球变化科学研究数据出版系统. 全球变化科学研究数据共享政策[OL]. https://doi.org/10.3974/dp.policy.2014.05 (2017年更新).
[12] 张兆明. 青藏高原Landsat系列卫星地表温度产品（2020）[Z]. 国家青藏高原科学数据中心, 2022. DOI: 10.11888/Terre.tpdc.272304.
[13] Wang, M. M., Zhang, Z. J., Hu, T., et al. A practical single-channel algorithm for land surface temperature retrieval: application to landsat series data [J]. Journal of Geophysical Research: Atmospheres, 2019, 124: 299–316.
[14] 程洁, 董胜越, 施建成. 中国区域1km无缝地表温度数据集（2002–2020）[Z]. 国家青藏高原科学数据中心, 2021. DOI: 10.11888/Meteoro.tpdc.271657.
[15] Xu, S., Cheng, J. A new land surface temperature fusion strategy based on cumulative distribution function matching and multiresolution Kalman filtering [J]. Remote Sensing of Environment, 2021, 254: 112256.
[16] Zhang, Q., Wang, N. L., Cheng, J., et al. A stepwise downscaling method for generating high-resolution land surface temperature from AMSR-E data [J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2020, 13: 5669–5681.
[17] Zhang, Q. Cheng, J. An empirical algorithm for retrieving land surface temperature from AMSR-E data considering the comprehensive effects of environmental variables [J]. Earth and Space Science, 2020, 7: e2019EA001006.